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LONG-TERM  GOALS 

The  development,  validation  and  application  of  robust  uncertainty  quantification  methods  to  ocean 
modeling  and  forecasting. 

OBJECTIVES 

This  project  explores  the  use  of  Polynomial  Chaos  (PC)  expansions  for  improving  our  understanding  of 
uncertainties  in  Ocean  General  Circulation  Models  (OGCM).  Given  adequate  initial  and  boundary 
conditions,  most  OGCMs  can  be  used  to  forecast  the  evolution  of  the  oceanic  state  consistent  with 
known  physical  laws.  Reliable  ocean  forecasts,  however,  require  an  objective,  practical  and  accurate 
methodology  to  assess  the  inherent  uncertainties  associated  with  the  model  and  data  used  to  produce 
these  forecasts.  OGCMs  uncertainties  stem  from  several  sources  that  include:  physical  approximation 
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of  the  equations  of  oceanic  motion;  discretization  and  modeling  errors;  an  incomplete  set  of  sparse  (and 
often  noisy)  observations  to  constrain  the  initial  and  boundary  conditions  of  the  model;  and 
uncertainties  in  surface  momentum  and  buoyancy  fluxes. 

Our  objective  is  the  development  of  an  uncertainty  quantification  methodology  that  is  efficient  in 
representing  the  solutions  dependence  on  the  stochastic  data,  that  is  robust  even  when  the  solution 
depends  discontinuously  on  the  stochastic  inputs,  that  can  handle  non-linear  processes,  that  propagates 
the  full  probability  density  functions  without  apriori  assumption  of  Gaussianity,  and  that  can  be  applied 
adaptively  to  probe  regions  of  steep  variations  and/or  bifurcation  in  a  high-dimensional  parametric 
space.  In  addition  we  are  interested  in  developing  utilities  for  decision  support  analysis;  specifically,  we 
plan  to  demonstrate  how  PC  representations  can  be  used  effectively  to  determine  the  non-linear 
sensitivity  of  the  solution  to  particular  components  of  the  random  data,  identify  dominant  contributors 
to  solution  uncertainty,  as  well  as  guide  and  prioritize  the  gathering  of  additional  data  through 
experiments  or  field  observations. 

APPROACH 

Our  approach  to  uncertainty  quantification  (UQ)  is  built  on  applying  PC  expansions  (Le  Maitre  and 
Knio,  2010)  to  investigate  uncertainties  in  simulating  the  oceanic  circulation  in  the  Gulf  of  Mexico.  We 
have  opted  to  use  the  HYbrid  Coordinate  Ocean  Model  (HYCOM)  as  our  simulation  engine  because  it 
has  been  developed  as  the  next  generation  model  for  the  US  Navy  and  has  been  adopted  by  NOAAs 
National  Center  for  Environmental  Prediction.  HYCOM  is  also  equipped  with  a  suite  of  sequential 
assimilation  schemes  that  will  be  used  to  investigate  how  UQ  may  be  beneficial  to  data  assimilation. 
Details  about  the  model,  its  validation,  and  sample  applications  can  be  found  in  Bleck  (2002);  Halliwell 
(2004);  Chassignet  et  al.  (2003,  2006)  as  by  visiting  http://www.hycom.org.  We  have  chosen  the  Gulf  of 
Mexico  for  our  experimentation  because  of  the  pre-existing  HYCOM  configurations  for  that  basin  and 
the  wealth  of  experience  accumulated  in  simulating  its  circulation,  the  availability  of  observations  for 
model  validation,  and  for  the  importance  of  that  region  to  the  social  and  economic  prosperity  of  the 
South  East  United  States.  In  the  following  subsections  we  present  the  main  ideas  of  the  PC  expansion 
before  we  summarize  our  efforts  to  date. 

PC  expansions  express  the  dependency  of  the  solution  on  the  uncertain  parameter  as  a  series  expansion 
of  the  form: 

p 

u(x,t,£)=Y,Mx,t)'i’k(£)-  (1) 

k= 0 

Here  u(x,  t,  £ )  is  a  model  solution  that  depends  on  space  x,  time  t  and  the  uncertain  parameters  £ ; 

*P is  a  suitably  chosen  orthogonal  basis;  and  u^xj)  are  the  expansion  coefficients.  In  our  approach 
to  UQ,  u  can  represent  a  variable  expressed  directly  in  the  model  such  as  sea  surface  temperature  or 
velocity  at  a  specified  point,  or  a  derived  quantity  such  as  the  mean  surface  cooling  under  a  hurricane 
track.  In  the  jargon  of  UQ  the  quantity  of  interest  u  is  often  referred  to  as  an  observable. 

The  choice  of  basis  function  is  dictated  primarily  by  the  probability  density  function  of  the  uncertain 
input  data,  /?(£),  which  enters  all  aspects  of  the  UQ  computations.  These  can  be  made  much  more 
efficiently  if  the  basis  vectors  are  orthonormal  with  respect  to  p(%).  Hence  the  basis  functions  are 
Legendre,  Hermite,  or  Laguerre  polynomials  when  the  input  uncertainty  is  described  by  uniform, 
Gaussian,  or  Gamma  distributions,  respectively. 
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The  computation  of  the  stochastic  modes  is  best  achieved  by  the  so-called  Non-Intrusive  Spectral 
Projection  (NISP)  method  since  we  would  like  to  avoid  modifying  the  original  OGCM  code.  Taking 
advantage  of  the  orthonormality  of  the  basis,  NISP  works  by  projecting  the  solution  u  on  the  basis 
function  via  inner  products,  and  by  replacing  the  integrals  with  quadrature  formula.  The  coefficients 
can  then  be  computed  simply  by  running  the  model  at  specified  values  of  the  uncertain  parameters, 
storing  the  desired  observable,  and  post-processing  via  a  simple  matrix-vector  multiplication.  No 
modification  to  the  OGCM  need  be  performed. 

The  investigative  team  at  Johns  Hopkins  University  consisted  of  Dr.  Omar  M.  Knio  (now  at  Duke 
University),  and  his  post-doctoral  associates,  Dr.  Alen  Alexanderian  and  Ihab  Sraj,  and  graduate 
student,  Justin  Winokur;  they  have  concentrated  on  advancing  the  technical  aspects  and  theoretical 
aspects  of  the  Uncertainty  Quantification  efforts.  Dr.  Mohamed  Iskandarani,  Ashwanth  Srinivasan  and 
William  C.  Thacker  (University  of  Miami)  have  focused  on  formulating  the  oceanographic  uncertainty 
problems,  the  modification  to  the  HYCOM  code  and  its  actual  execution,  and  the  preparation  of  the 
necessary  data  to  carry  out  the  research  agenda.  Dr.  Matthieu  Le  Henaff  assisted  us  with  the  Gulf  of 
Mexico  configuration,  and  we  have  held  discussions  with  Dr.  Francois  Counillon  as  to  the  applicability 
of  PCs  in  Ensemble  Kalman  Filters  based  data  assimilation. 

WORK  COMPLETED 

Our  work  proceeded  along  several  parallel  tracks:  studying  the  uncertainty  stemming  from  nesting 
boundary  conditions,  investigating  parametric  uncertainty  in  HYCOM  mixed-layer  model  and 
wind-drag  bulk  formula,  assessment  of  uncertainties  in  an  oil  spill  model,  and  developing  an  efficient 
procedure  to  handle  stochastic  dynamical  systems  with  multiple  time  scales. 

Our  initial  experimentation  with  UQ  revolved  around  the  impact  of  nesting  boundary  conditions 
uncertainty  on  the  ocean  circulation  forecasts  in  the  Gulf  of  Mexico.  In  addition  to  the  development  of 
the  non-intrusive  HYCOM  and  UQ  machinery  to  carry  out  this  task,  we  proposed  an  original  approach 
to  cast  the  problem  within  a  PC-framework  so  that  only  two  stochastic  parameters  are  needed.  The 
nesting  data  was  decomposed  into  space-time  Empirical  Orthogonal  Functions  and  the  two  dominant 
modes  were  perturbed.  An  ensemble  of  49  simulations  were  needed  to  characterize  the  entire 
uncertainty  space.  Figure  1  shows  the  evolution  of  the  17  cm  surface  height  field  for  these  49  runs,  and 
figure  2  shows  PC-estimates  of  the  standard  deviation  at  various  times.  An  article  describing  this  work 
has  been  submitted  for  publication  (Thacker  et  al.,  201 1). 

The  Deep  Water  Horizon  oil  spill  afforded  us  the  opportunity  to  use  UQ  for  a  pressing  environmental 
problem,  and  to  help  understand  the  fate  and  potential  impact  of  the  discharged  oil.  The  simulation 
combined  two  nested  HYCOM  models  (at  2  km  and  900  m  resolutions)  and  a  3D  oil  spill  model  that 
included  the  effects  of  advections  by  ocean  currents  and  waves,  and  the  effects  of  various  physical, 
chemical  and  biological  processes.  The  uncertainty  in  the  initial  distribution  of  oil  droplet  size  and 
compositions  were  propagated  through  the  model  using  the  UQ  machinery  and  the  results  are  illustrated 
in  figure  3.  It  should  be  noted  that  this  work  required  a  herculian  effort  to  marshal  and  synchronize  a  set 
of  heterogeneous  and  distributed  high  performance  computing  resources  to  carry  out  the  simulations; 
this  accomplishment  has  been  recognized  in  a  talk  and  a  publication  in  a  prestigious  journal  Srinivasan 
etal.  (2010). 
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The  parametric  study  focused  on  quantifying  ocean  model  uncertainties,  associated  with  mixed  layer 
and  air-sea  momentum  exchange  parametrizations,  when  the  ocean  is  forced  by  hurricane-strength 
winds.  The  wind-stress  and  the  mixed  layer  parameterization  modulate  the  ocean  Sea  Surface 
Temperature  (SST)  response  to  hurricane  forcing,  a  key  parameter  in  determining  the  exchanges  of 
energy  and  momentum  between  the  ocean  and  atmosphere,  and  in  impacting  hurricane  intensity 
forecast.  The  setting  chosen  for  our  experiments  is  the  circulation  in  the  Gulf  of  Mexico  during  the 
passage  of  hurricane  Ivan  from  Sep  9-16  2004.  We  used  PC  expansions  to  systematically  characterize 
the  entire  response  surface  of  the  ocean  model  to  parametric  uncertainties  in  its  subgrid 
parameterizations.  Four  parameters  were  perturbed  in  this  study:  three  KPP  parameters  (The  critical 
Richardson  number  controlling  the  definition  of  the  ocean  mixed  layer,  and  the  background  viscosity 
and  diffusion),  along  with  the  wind-drag  coefficient.  A  total  of  385  simulations  were  necessary  to 
explore  the  entire  uncertainty  space  and  determine  all  the  stochastic  modes;  various  consistency  checks 
were  done  to  confirm  that  the  PC  representation  remained  accurate.  Figures  4-6  illustrate  the  statistical 
output  that  can  be  mined  from  these  simulations.  Figure  4  shows  the  mean  mixed-layer  depth  as 
estimated  by  the  PC  expansions.  Figure  5  shows  the  evolution  of  the  SST  probability  density  function  at 
a  NOAA  buoy  (located  at  85.1  W  and  26.1  N);  the  latter  shifts  to  cooler  temperature,  as  expected  in  the 
hurricane  wake,  and  broadens  with  time.  Figure  6  illustrates  how  sensitivity  analysis  metrics  can  be 
used  to  identify  the  main  contributors  to  SST  uncertainty  throughout  the  Gulf  of  Mexico. 

A  preconditioning  approach  was  developed  to  improve  the  efficiency  of  representing  uncertainty  in 
dynamical  systems  with  PC  expansions.  The  approach  is  based  on  the  definition  of  an  appropriate 
multiscale  stretching  of  the  individual  components  of  the  dynamical  system,  and  which  can  robustly 
recovery  the  unsealed  transient  dynamics.  The  stochastic  modes  are  then  obtained  through 
non-intrusive  spectral  projections  of  the  stretched  measures.  Implementation  of  the  present  approach 
was  illustrated  through  application  to  a  chemical  system  with  large  uncertainties  in  the  reaction  rate 
constants.  Computational  experiments  showed  that,  despite  the  large  stochastic  variability  of  the 
stochastic  solution,  the  resulting  dynamics  can  be  efficiently  represented  using  sparse  low-order  PC 
expansions  of  the  stochastic  multiscale  preconditioner  and  of  stretched  variables.  We  will  explore 
extending  these  preconditioning  ideas  to  oceanic  simulations  in  the  near  future. 

The  work  associated  with  this  project  has  been  publicized  at  several  conferences,  workshops,  seminars 
and  invited  talks,  including: 

•  the  workshop  on  Uncertainty  Quantification  for  Multiscale  Systems,  Baltimore  MD,  July  20  2010; 

•  the  workshop  on  Mathematical  Theory  and  Modelling  in  Atmosphere-Ocean-Science,  Oberwolfach 
Germany  08/08/2010-08/14/2010; 

•  “Many  Task  Computing  for  Modeling  the  Fate  of  Oil  Discharged  from  the  Deep  Water  Horizon  Well 
Blowout”,  MTAGS  2010,  3rd  IEEE  Workshop  on  Many-Task  Computing  on  Grids  and 
Supercomputers,  New  Orleans  LA,  11/15/2010; 

•  the  Layered  Ocean  Model  workshop,  Miami  FL,  02/07/201 1-02/09/2011; 

•  the  SIAM  Conference  on  Computational  Science  and  Engineering,  Reno  NV, 
02/28/2011-03/04/2011; 

•  the  SIAM  Conference  on  Mathematical  and  Computational  Issues  in  the  Geosciences,  Long  Beach 
CA,  03/21/201 1-03/24/2011. 

•  “Probabilistic  Methods  for  Uncertainty  Quantification:  Recent  Developments  and  Implications  to 
Extreme  Scale  Computing,”  DOE-SCGF  Annual  Meeting,  July  18,  2011. 

•  “Stochastic  Multiscale  Approaches  in  Fluid  and  Thermal  Systems,  KAUST,  May  22,  2011. 

•  “Stochastic  Multiscale  Approaches  in  Fluid  and  Thermal  Systems,”  Masdar  Institute,  May  8,  2011. 

•  “Multiscale  Simulations,  Uncertainty  Quantification  and  Data  Challenges,”  University  of  Lugano, 
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March  28,2011. 

•  “Multiscale  Paradigms  in  Reactive  Solids  and  in  Fluids,”  Duke  University,  March  22,  2011. 

•  “Spectral  Approaches  of  Uncertainty  Quantification  in  Complex  Dynamical  Systems,”  UC  Berkeley, 
February  16,  2011. 

•  “Spectral  Approach  to  Uncertainty  Quantification  in  Complex  Dynamical  Systems,”  plenary  lecture, 
2010  MetStroem  Program  Meeting,  FU  Berlin,  October  28,  2010. 

RESULTS 

The  Empirical  Orthogonal  Function  decomposition  proved  to  be  an  effective  method  to  recast  the 
uncertainty  quantification  problem  into  a  PC-based  practical  framework.  A  similar  decomposition  can 
be  envisioned  for  perturbations  to  initial  conditions.  This  is  an  on-going  investigation. 

A  number  of  scripts  have  been  developed  to  coordinate  the  launching  of  the  HYCOM  ensemble  and  the 
collection  of  the  output  for  UQ  post-processing.  In  addition,  Smolyak  nested  integration  rules  were 
coded  to  mitigate  the  growth  in  the  size  of  the  ensemble  as  the  number  of  stochastic  variables  increases. 
This  nested  quadrature  is  a  good  basis  upon  which  to  build  adaptive  sampling  techniques. 

The  PC-premise  of  casting  the  probabilistic  uncertainty  as  an  approximation/interpolation  problem  is 
paying  off  in  a  couple  of  ways.  First  it  allows  for  self-consistency  checks  as  to  the  error  committed  in 
the  approximation,  and  second  it  allows  the  expansions  to  act  as  a  cheap  emulator  for  the  model.  Thus, 
histograms  of  output  uncertainties  can  be  built  cheaply.  More  importantly  various  variance  and 
sensitivity  analysis  metrics  can  be  devised  to  assess  which  of  the  uncertain  input  parameters  have  the 
greatest  impact  on  the  output  uncertainties.  This  is  illustrated  in  figure  6  which  shows  how  the  fraction 
of  variance  in  SST  that  can  be  attributed  to  the  uncertain  input  parameters;  here  the  background 
diffuvisity  is  the  dominant  factor  in  the  absence  of  strong  wind  forcing  whereas  the  drag  coefficient 
takes  over  in  the  hurricane  track  region. 

IMPACT/APPLICATIONS 

The  present  project  presents  an  approach  to  characterize  the  entire  response  surface  of  an  ocean  model 
to  uncertainties  in  its  input  data.  This  has  implications  for  the  field  of  data  assimilation,  particularly  for 
ensemble  Kalman  filter  based  approaches.  The  methodology  developed  here  will  be  of  use  either  for  the 
efficient  update  of  the  covariance  matrices  and/or  quantifying  the  errors  incurred  by  small  size 
ensembles.  We  plan  to  explore  these  ideas  in  the  near  future. 

TRANSITIONS 
RELATED  PROJECTS 

Dr.  Ashwanth  Srinivasan  was  partially  supported  by  an  NSF-RAPID  grant  (NSF  OCE- 1048697)  for  his 
work  on  the  oil-fate  model. 
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Figure  1:  Sea  Surface  Height  contours  from  the  ensemble  HYCOM  runs  at  the  quadrature  points. 
The  panels  show  the  state  of  the  17  cm  contour,  used  as  a  proxy  for  the  Loop  Current,  at  different 
times.  The  contour  from  each  deterministic  run  of  HYCOM  (realization)  is  plotted  using  a  unique 

color 
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Figure  2:  Standard  Deviation  of  the  Surface  Height  field  obtained  from  the  Polynomial  Chaos 
Expansions.  The  panels  show  the  evolution  of  the  uncertainty  every  150  days.  The  standard 
deviation  is  computed  using  a  6  order  polynomial  expansion. 
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Figure  3:  The  upper  figure  shows  the  surface  slick  simulated  by  the  oil  model  for  June  11,  and 
compares  well  with  satellite  based  composites.  The  lower  figure  shows  the  subsurface  fraction  of  oil 
between  1000-1400  m  as  a  function  of  time.  The  error  bars  are  based  on  input  uncertainties  on  the 
oil  droplet  size  and  the  oil  composition.  The  results  suggest  that  about  80  %  of  the  oil  has  been 
removed  from  the  water  column  by  evaporation  and  degradation  and  the  uncertainty  analysis 

assigned  a  20%  uncertainty  to  this  estimate. 


10 


Expocted  value  ol  mid  at  t »  150 


Expected  value  ol  mW  at  t  =  162 


© 

3 


Longitude 


-95  -90 

Longitude 


Figure  4:  Time  evolution  of  the  expected  Mixed  Layer  Depth.  The  red  curve  and  the  black  dot  refer 

to  Ivan  ’.s'  obeserved  track  and  to  its  center,  respectively. 


Figure  5:  Distribution  of  SST  at  the  buoy  location  over  time  ( generated  by  sampling  the  PC 
expansions  with  sample  size  of  106.)  Note  that  as  the  hurricane  nears  the  Gulf  (t  «  90),  enters  the 
Gulf  (t  ~  120)  and  gets  close  to  the  buoy  (t  ~  150)  we  see  the  cooling  of  water  through  the  shift  of  the 
distribution  to  left.  Also,  note  the  emergence  of  long  tails  for  the  distributions  stretching  toward 

cooler  temperatures. 
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Figure  6:  SST  sensitivity  analysis  with  Hurricane  Ivan,  before  (t=90  hr)  and  after  (t=147  hr)  the 
hurricane  enters  the  Gulf  Here  we  look  at  T$  (  sensitivity  due  to  background  diffusivity)  and  T4 
( sensitivity  due  to  stochastic  wind  drag  coefficient).  Background  diffusivity  accounts  for  most  of  the 
variance  in  the  absence  of  the  hurricane  while  the  drag  coefficients  has  the  most  influence  in  the 

region  of  strong  wind  forcing  ( right  column). 
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